Developing a computational model of blood platelets 
with fluid dynamics applications 

Vijay Viswanathan 
March 2012 

vi j ayv@andrew . emu . edu. 



1 



Table of Contents 

Contents 

1 Introduction 

2 Methods 

3 Results 

4 Analysis 

5 Conclusion and Future Work 



Vijay Viswanathan 




Figure 1: As a graphical aid, here is an image [T], not taken by the author, which demon- 
strates platelet aggregation around a red blood cell. 

1 Introduction 

Note: this paper was originally submitted to the Siemens Competition for Math, Science, 
and Technology, so all conventions follow Siemens guidelines. 

This paper considers studying blood platelets using mathematical analysis and meth- 
ods. Platelets, also known as thrombocytes, play a key role in blood clotting. The role 
of these entities in strokes, myocardial infarctions, and coronary artery disease add to the 
importance of blood platelets. The chemical explanations of blood platelet activation and 
coagulation are rather well understood. Still, the mechanics behind platelet aggregation is 
largely unknown. According to researcher Constantine Pozrikidis of the University of Mas- 
sachusetts at Amherst, "While the biochemical origin of the adhesion kinetics has been well 
characterized, the non-spherical platelet shape has discouraged the mathematical modeling 
of the adhesion process by elementary methods of particulate hydrodynamics." [2] The 
fact is that the shape of platelets is derailing the hydrodynamic simulation efforts of blood 
platelet activation, aggregation, and adhesion, and this is hampering scientific progress in 
hematology. 

Due to the microscopic nature of blood platelets (the mean platelet volume for an av- 
erage 65-year old is only 10 femtoliters [3] - 100 millionths of a cubic millimeter) accurate 
laboratory results involving thrombocytes are very difficult to achieve. A newly developing 
method to approximate these results is computational fluid dynamics (CFD). Using meth- 
ods of CFD, complex experiments involving fluid flows, turbulence, and interactions between 
liquids, gases, and solids can be solved using numerical algorithms. However, as Pozrikidis 
states, mathematical and hydrodynamic models of platelets are generally undeveloped, and 
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this hinders the growth of CFD methods with regards to thrombocytes. 

When blood platelets enter the coagulation cascade and aggregate with other platelets, 
they undergo the process of platelet activation by thrombin. Morphologically, this makes 
the shape of the platelet shift from an ellipsoid-like object Q to a convex, irregular body. 
While unactivated platelets are easily representable and have been modeled in the past, 
activated platelets are more challenging to computationally work with. The reasons for this 
are twofold: activated platelets have a highly irregular shape, and due to the nature of the 
platelet aggregation process, very few mechanical results of individual platelets have been 
found. 

In developing a CFD simulation, either the geometry of the solids and fluids in question 
or more intrinsic characteristics (e.g. viscosity of fluids, density of solids, etc.) must be 
known. Owing to the complex shape of platelets, the former path becomes computationally 
inefficient. Similarly, very few existing CFD simulators offer support for shapes which lack 
a simple functional representation. 

Regarding the intrinsic characteristics of the objects in a simulation of platelets, one pri- 
mary unknown variable is the drag force caused by the platelet. Again, the blood platelet 's 
odd shape makes evaluating the drag coefficient (c<j) rather difficult from a theoretical stand- 
point. 

In response to the gaps present in both pathways towards a CFD simulation involving 
thrombocytes, I decided to answer the question of how one can develop a computational 
model of blood platelets for fluid dynamics purposes. 

This process took three stages. The first step in the process of developing the model 
was outlining the mathematics behind the platelet geometry, both in two-dimensions and 
in three-dimensions. Next, I developed a method to computationally "match" any platelet 
image to the analytical expressions that represent the shape of the platelet. Finally, the 
platelet images were utilized to attain the drag coefficient for each blood platelet analyzed. 
This could become extremely useful if scientists hope to design personalized devices (like 
stents) or if pharmacists hope to create computationally designed antiplatelet drugs for 
patients based on their specific platelet characteristics. This is currently unachievable, 
but up-and-coming simulations using the knowledge enumerated in this paper could help 
this dream. Hopefully accurate CFD simulations could feasibly replace arduous laboratory 
experimentation. 
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Figure 2: Here is the first approach of the two-dimensional fitting phase, where for every 
section of the edge of the image, polynomial interpolation created a polynomial approxima- 
tion. 



Note: all programming was performed in C++, primarily using default libraries in Visual 
C++ 2010 Express, as well as the CML T Graphics Library |5j. Any functions which do 
not belong to default libraries or CMU graphics were written by the author of this paper. 
Similarly, all two-dimensional graphing was performed on author-generated software while 
all three-dimensional graphing was done on Mathematica. 
The methodology of this project lies in three phases: 

Phase 1: Two-dimensional modeling. During the process of achieving a workable 
three-dimensional model, the first step (for simplicity) was modeling in two-dimensions. 
Initially, I experimented with the feasibility of taking a platelet image and interpolating 
the points on the edges of the image with low-order polynomials. However, it became clear 
that while these formulations were rather faithful to the initial platelet image they were 
both unnecessary and inefficient. For one, there existed some portions of the shapes (such 
as vertical lines) which could not be modeled with polynomials with the same variables. 
Additionally, such complex curve-fitting could not be generalized into an effective framework 
with a set number of parameters, and the time effectiveness of precise interpolation in three 
dimensions would seriously stifle the relevance of such an algorithm. 

Next, I developed a general framework using the classic curve in polar coordinates known 
as the "polar rose" but generalized it to fit the typical model of a blood platelet. A platelet in 
a simplified two-dimensional visual looks like an invisible circle with projections of different 
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Figure 3: Above is the second approach in the two-dimensional phase, where using the 
generalized polar rose model, a platelet-like object was generated 

lengths and widths emerging radially from it. To make the polar rose have this inside circle, 
a "generalized polar rose" which randomly varies the parameters of individual rose petals 
was created. Furthermore, due to the unavailability of existing software capable of graphing 
piecewise, polar functions, I designed and implemented my own graphing utility capable of 
randomizing parameters and graphing piecewise-defined polar roses. 

Phase 2: Three-dimensional expressions. The initial stage of this phase was to 
implement the equation of an ellipsoid. This is the core of the platelet model, and projections 
emerge outwards from the surface of this ellipsoid. 

The second stage for the three-dimensional mathematical framework of the platelet mor- 
phology was developing the equations of the projections. To improve smoothness between 
these surface projections and the underlying discoid structure, a logarithmic function is 
used for the surface projection. By solving a simple differential equation, the function rep- 
resenting the projections can be made to fit with the ellipse in a manner which maintains 
differentiability for the most part. In order to generate a three-dimensional surface from 
this 2D curve, this function was revolved around the z-axis and surfaces of revolution were 
employed. Additionally, parameters that allow the changing of the length and width of each 
projection were exercised. 

In order to allow these projections to emerge in all directions rather than exclusively 
upwards, different approaches were taken. Initially, spherical coordinates were used as a 
simply mode of rotating the upwards projection to many different direction. However, after 
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Figure 4: Pictured is a sample model of a platelet using the framework developed in this 
paper. 

converting the z(x, y) function to an r{9,(j)) function in spherical coordinates, the function 
could no longer be solved explicitly (i.e. it was not possible to isolate the r variable from 
the other two variables). By consequence, while rotations would have been easy to perform, 
numerical algorithms would be required to evaluate the function in spherical coordinates. 
Then, cylindrical coordinates were engaged to avoid the complications found in spherical 
coordinates. However, once the explicit form was achieved in cylindrical coordinates, rota- 
tions could only be done along one angle, and this made cylindrical coordinates impractical. 
To maintain simple and explicit functions, ultimately rotational matrices had to be applied. 
This rotational matrix was a reduced version of the traditional three-angle rotation matrix 
in three dimensions, because one of the Euler angles (the "roll" angle) was not relevant to 
the rotation of projections. 

Phase 3: Computer program and drag coefficient calculation 

Note: Specific figures and screenshots referred to are in the Illustration section. 

This portion of the work required writing a computer program in CH — h which facilitated 
the process of matching any platelet to a mathematical representation. The platelets were 
inputted with six different, orthogonal views. For example, an object in three-dimensional 
Cartesian space with variables (x,y,z) could be understood by viewing the images through 
the (x,y) plane, through the (x,z) plane, through the (y,z) plane, and through the view 
from the opposite face of each plane. Additionally, due to the difficulty and inaccuracy of 
completely automating the process, data by the user was manipulated into a model. 

Furthermore, the specific dimensions of the ellipsoid resting within the projections had 
to be found. 

The last endeavor of this phase was to find the drag force created against fluid flow 
by blood platelets. In this case, there was no existing literature describing the drag force 
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exerted by an isolated platelet in a flow, thus necessitating the use of a theoretical basis for 
drag coefficient calculations. Tran-Cong, Gay, and Michaelides [B] in 2003 created a method 
to approximate the drag coefficient (c<j) for any object in a fluid flow (with the condition 
that the flow's Reynolds number had to be between 0.f5 and 1500. 

Most of the variables in the Cd formula were applicable to all platelets and could be 
found in recent literature. My written code used algorithms to determine two of the values 
required in this theoretical Cd formula. 

3 Results 

Results from Phase 1: Fig. 2 displays a picture of a platelet image with the polynomial 
interpolations graphed over it. This shows the primary result for this stage, when the goal 
was to curve-fit every section of the image's border to a polynomial. 

The next result gained was the two-dimensional model using the idea of a generalized 
polar rose. The equation for the generalized polar rose is r = a ■ cos(b9 + 4>) + c, where a 
indicates projection length, b indicates projection width, <j> rotates the curve, and is the 
variable for the function. The written program randomized the lengths and widths while 
connecting the different projections in a smooth manner. This was done by first splitting 
the total, 360° curve into subsections of different widths (done by reducing the b value in 
the generalized polar rose equation so that the quantity of petals in the defined interval is a 
positive integer). Then, within each section of identical width, individual petals were given 
random heights. This simple method created a random platelet each time the program was 
run. Figs. 3 and 4 display sample random drawings made by the program to approximate 
the two-dimensional appearance of a platelet. In order to display these piecewise polar 
functions, a simple polar grapher with support for separate pieces on different intervals 
was developed. In Figs. 3 and 4, the randomized equations are shown as graphed on this 
piecewise, polar-curve graphing utility. 

Results from Phase 2: The initial stage of this phase was to implement the equation 
of an ellipsoid. The simple, Cartesian equation for an ellipse with semidiameters a, b, and c 

9 9 9 

x y z 

is — ^ + H — 77 = 1. This is the core of the platelet model, and projections emerge outwards 
a 2 b z cr 

from the surface of this ellipsoid. 

In developing the three-dimensional projection functions, a two-dimensional basis func- 
tion had to be made. The chosen function type was the logarithmic function. This decision 
was determined after solving the differential equation /'(0) = 1/0 (made so that when ro- 
tated the point at wouldn't be a cusp, but rather a smooth end). One possibility for 
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Figure 6: Attached is another generalized, randomized polar rose, of the equation r = 



f'(x) was the function f'(x) = 1/x. Integrating on both sides, the resulting function is 
/ f'(x) dx = J 1/x dx => f(x) = log(x) + C. In order to generate a three-dimensional sur- 
face from this 2D curve, this function was revolved around the z-axis. As a method finding 
the equation for the surface of revolution, the radius of the circular level sets of the surface 
was set to vary with y in a logarithmic manner. This level set method served as the process 
to find the projection's surface equation. Furthermore, to adjust the length (denoted by I) 

and width (denoted by w) of the projections and for the projections to be outwardly project- 

(log[—x + I}) 2 

ing, the specific two- variable function chosen was y 2 + z 2 = — ■ . The projection 



(log[—x + 1}) 2 

surface that resulted can be represented by the original function z = ±\/ y 2 . 

V w 

For these projections to be fitted to face any direction, rotational matrices had to be 
applied. To determine the general form for rotating the function by two angles, 6 and <$> 
(colloquially known as pitch and yaw in terms of Euler angles) the rotation matrix must be 
multiplied with the vector of the original parametric function. The matrix shown below is 
the rotation matrix when ip, the third Euler angle representing the "roll" motion, was set 
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to 0. 

cos9cosip simp sinOcosip 

cosQsimp cosip sinOsimp 

—sinQ cos9 

By consequence, the matrix multiplication between the rotation matrix and the vector 
representing the function is: 



cosQ sincpsinO cos(psin0 

cos<p —sirup 

—sinO —sincpcosO cos(pcos6 



z{x, 9, <p, I, m)- 



x 

y 



(log[-x + l}Y 
w 



ucosQ + vsincpsinO + cos(psin9 



vcos0 — sine 



(log[-x + l)) 2 



-xsinO — ysincpcosO + coscpcosO 



{log[-x + l]f 



w 



y 



z{x 1 9 1 (p,l,m) represents a projection with length I, width m, and rotated angles 9 and 
<p from the z-axis. 

Results from Phase 3: Note: Specific figures and screenshots referred to are in the 
Illustration section. 

The next step was developing the computational aspects of the program. Specifically, the 
program was supposed to fit the developed mathematical model and match the parameters 
to any inputted platelet image (represented with six different two-dimensional views of the 
same three-dimensional object). The program had to gather data on the length and width 
of each platelet, the dimensions of the underlying ellipsoid beneath the projections, the area 
of the platelet from one view, and the perimeter of the platelet from the same view. 

The application of the six-view method for analyzing the platelet image became difficult 
as these specific views have never been taken before. As a result, I approximated the shape of 
a platelet by taking the existing mathematical model and randomizing the parameters. This 
created a varied shape which was a suitable replacement for the blood platelet morphology. 

The first action was developing an interface on which a user could aid the computer in 
the detection of the projections. Automatic detection of these projections took too much 
runtime and was not sufficiently accurate given the typical low-contrast of the images, so 
significant user input was required. Initially, the six sides of a blood platelet were juxtaposed 
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Figure 9: Comparison between the initial platelet image and what was later detected. Issues 
with a third-party software account for the apparent difference, but fundamentally both 
images have the same elements. Also, the projections are not correctly bounded in the 
second image, as based on the mathematics of the model. That is, projection A should be 
only show the right half side, projection should only show its left half side, etc. 

on the same graphics window. Then, the user clicked on the point representing the endpoint 
of the same projection in three different views - this way x, y, and z coordinates could be 
shown for each endpoint. In the same function, a utility for determining the width of each 
platelet was defined by allowing the user to click on the endpoints of each projection's base. 

Moreover, the dimensions of the underlying ellipsoid had to be discovered. Again, manual 
detection with automatic computation was favored, and the user clicked on the endpoints 
of each of the three semidiameters of the ellipsoid. As a result, the a, &, and c components 

2 2 2 

x y z 

in the ellipsoid equation, + H — ^ = 1, were identified. In order to compare these 

a 1 tr c z 

equations with real-life data, it was found that the length of 1 pixel in the graphics window 
corresponded to 10.62 nanometers, and this leads to accurate predictions of platelet size. 

The final method in this work was the determination of the drag coefficient for any 
inputted blood platelet. In 2003, Tran-Cong, Gay, and Michaelides j6| developed a formula 
describing the drag coefficient of an irregular object. The authors also mentioned that this 
formula is consistent with experimental results for 0.15 < Re < 1500. Their formula is 
shown below in Equation 1. 



24 d A r 0.15 eU p xo.6871 
Re d n Jc a n 



0.42(^) 2 



^[1 + 4.25 x 10 4 ( C ^i?e)- 1 - 16 ] 



(1) 



In Equation 1, Re is the Reynolds number of the flow, and for blood flow in a human 
artery, the average Reynolds number is 250. d n is the nominal diameter of the object, and 
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/ 6V 7 " 

d n = (\ — ) 3 > where V is the particle volume The specific value for V that we use is 10.0 

pLA~ 

femtoliters [3]. d a is the volume-equivalent-sphere-diameter. d a — y — - (where A p is the 

projected area of the object in the direction of motion). Lastly, c, the particle circularity, 

■nda 
K 

in the direction of motion. 



can be found with the equation c = —f^-, where P p is the projected perimeter of the object 



The program that was developed found both A p and P p , and the remaining parameters 
either depended on A p and P p or were found in other publications. The drag coefficient 
calculated for the platelet shown in Fig. 4 and 5 was 2.874. This corresponds with drag 
coefficients for more hydrodynamic shapes, thus supporting the claim. 

Finally, to make these results applicable to real medical technology, the scale between 
the length of one pixel in the graphics component and physical length of platelets had to 
be found. Using the A p variable, sqrtA p 3 roughly represents the volume of an object with 

3 

all sides of area A p . Comparing WA P with the mean platelet volume (assumed to be 10.0 
femtoliters), are roughly 0.0000043 femtoliters per voxel. These means4.28 x 10 -21 liters 
per voxel. By taking the cubic root of both 4.28 x 10 -21 and the liters per voxel unit, the 
result is that the length of 1 pixel roughly corresponds to 1.62 x 10~ 8 . By this measure, the 
average length of 1 pixel of the platelet is roughly 16 nanometers. 



4 Analysis 

Going into this project, I had expected to design a basic mathematical model which would 
then go into a coarse-grained molecular dynamics model of platelet activation. However, 
my model was not discrete enough for a dissipative particle dynamics approach to coarse- 
grained molecular dynamics (MD). Moreover, were my project to become a component of 
an actual MD simulation, the runtime of the simulation would literally take years, even if 
the processes were parallelized on high-performance computing clusters. 

The transition from a computational model with MD applications to one with primarily 
CFD results was rather smooth given the nature of my model. The equations that I provided 
were continuous and, for the most part, smooth wherever defined. While my use of the 
logarithmic function in a model that hoped to be differentiable could be somewhat alarming, 
log(x) is differentiable on all orders for all values on which log(x) is defined. By consequence, 
the surface of the geometric platelet model is almost completely differentiable. 

My findings could not be referenced with external data because of a lack of published 
material on the fluid mechanics of thrombocytes. While some sources do reference the drag 
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force of platelets in their calculations [7J, never do they spefically mention values or equations 
for calculating the drag forces and coefficients. Similarly, Karniadakis et. al.[S] mention drag 
coefficients of platelets as generally varying inversely with the Reynolds number, but there 
is an absence clear documentation of Cd values. This chasm in knowledge is clear when 
looking at some influential blood platelet simulations. For example, Eckstein and Belcagem 
[Hj describe a model of platelet transport with drift and diffusion, but they also write, "No 
extra drag or interaction with the wall (beyond that included as a part of the drift) was 
included," with respect to platelets in their simulation. 

Rather than being supported by these other findings, my work enhances these findings. 
For example, in the poster by Karniadakis et. al. [8], more rigorous calculations of the drag 
coefficient could be computed without significant experimental issues. Similarly, in Eckstein 
and Belcagem's paper, they could improve the scope of their work by including not only 
drift and diffusion in the model, but also drag. Many avenues of research can be expanded 
upon with this computational method of discovering drag, because this reduces and in the 
case of blood platelets eliminates the need for many experimental procedures relating to 
these calculations. 

The expected result was for the mathematical model be reasonably faithful in its fit 
to a real platelet while not impeding on the model's efficiency. However, due to the lack 
of sufficient platelet images, this could not be tested - only the solution could be given. 
Still, when comparing the approximated platelet (generated by randomizing an ellipsoid 
with projections around it), while initially the result looks much busier than the original 
input, closer analysis shows that this effect is due to the different scales required for each 
parameter. The parametric equations have two parameters, u, and v, but in Mathematica 
the two parameters are always treated in the same interval for all parametric surfaces in 
the same graphing window. As a result, this does not account for the different intervals of 
graphing required for each surface (u and v should only be defined for each projection until 
the very tip of the projection can be seen - after that surface extraneities emerge). 

When varying different approximations of platelets as the input of the program, there 
was generally strong consistency. Only if the variation between projections was very strong 
would real issues arise, and even that was largely due to issues with the Mathematica 
graphing windows rather than anything else. The equations themselves were correct based 
on comparison between the resultant equations and the input equations that defined the 
"platelet" which served as the input for my code. 

Ultimately, this work is unique and is a contribution to science because it uses easy-to- 
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use computational methods to build on very difficult experimental problems. For example, 
stents to prevent infarctions could be designed using CFD software such that for a person 
with a rare blood disorder, the safety and efficiency of the stent could be assured. Similarly, 
my mathematical model could be used to develop normative criterion for "safe" platelet 
structure. Patients not falling completely into these normative intervals could be investi- 
gated further and diagnosed. For example, patients whose average projection length exceeds 
the national average would require further examination. 

5 Conclusion and Future Work 

My work accomplished the two goals that were set out. I developed analytic expressions 
representative of blood platelets in varying stages of activation and developed code that 
fit these expressions to faithfully represent any given blood platelet. As a result of this 
work, future computational fluid dynamics simulations involving blood platelets in laminar 
fluid flows as well as simple turbulent flow. Moreover, because this work provides two 
different bases for improved simulations (both a geometric model and a solution for the 
drag coefficient), it provides a solution for simulations on two platforms. 

My conclusions are supported by the results in the report alone because of the chasm 
between what results are computationally possible and what results are achievable through 
actual experimentation. Experimentally, the small size of blood platelets and the rapid 
changes that occur to thrombocytes hinder efforts to calculate platelet mechanics. As such, 
specific results, like the average drag coefficient of an activated platelet, are not available 
to support my current claims. The only method that could somewhat verify my drag 
calculations is by comparing my results to results for other, distinct shapes. For example, if 
the drag coefficient for a triangular prism pointing away from the direction of flow was 1.14 
[TU] then the drag coefficient for a blood platelet would have to be substantially greater due 
to the highly non-aerodynamic or hydrodynamic shape of the platelet. 

However, my methods are certainly not infallible. Most of the testing for my program 
relied on randomly-generated, analytical models. Therefore, my program fitted a model to 
an existing image of a model. Still, because of the variety of blood platelet shapes, this should 
not have caused major problems. Furthermore, some of my algorithms (such as my algorithm 
to find the perimeter of irregularly shaped images) could have been overestimations or 
undcrcalculations. 

Experiments that would have to be performed to refine the calculations would necessitate 
the experimental testing of my theoretical results. When molecular biological techniques, 
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microfluidics, and techniques of microscopy become sufficiently refined, accurate fluid dy- 
namic calculations of blood platelets in blood flow will become achievable. The theoretical 
results in this paper could then be checked with the experimental results and statistical 
measures could potentially ensure significance. Likewise, in the future other mathematical 
frameworks to fit platelets and other related fluid dynamics calculations could be performed 
and referenced with the current theoretical results to check for accuracy. 

Additionally, to improve this research, the detection of the characteristics of each platelet 
(like the lengths and widths of projections and the dimensions of the underlying ellipsoid) 
must be efficiently automated. 

This research aimed to be more valid than any other currently published, computational 
results. In fact this seems to be the case, as an exhaustive literature review found no clear 
mathematical model for blood platelet morphology and results for the drag coefficient of 
blood platelets are nonexistent. 

My future for this work is that simulations can be used to develop the way medicine is 
approached. Through this paper I attempted to patch up many of the holes that inhibit 
accurate and efficient computational fluid dynamics simulations. In the future there would 
be two primary concepts which could be expanded upon. Firstly, more sophisticated com- 
puter vision techniques could be applied in the creation of a model for medical images. This 
would completely automate the process and improve general applicability. Second, similar 
methods of morphological model matching can be used for other aspects of medicine and 
biology. 
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